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ABSTRACT: 



This report presents an operational computer program for the 
minimal-time routing of single ships. Although written specif- 
ically for VC2AP3 and VC2AP2 vessels operating in a described 
area of the north Pacific ocean, the program can be modified 
easily to provide routes for other type vessels in any ocean 
area of the northern hemisphere. The method of incorporating 
weather forecasts into the preparation of minimal-time ship 
routes is described, and possible future developments are dis- 
cussed. 
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I . Introduction . 

Section II, by Prof. G. J. Haltiner, describes the scheme 
of incorporating weather forecasts into the preparation of mini- 
mal time ship routes used in this report, and discusses possible 
future developments. The remaining Sections of the report extend 
the previous work of Profs. W. E. Bleick and F. D. Faulkner, [l] 
and [ 2 ], by describing an operational computer program for the 
minimal-time routing of VC2AP3 vessels in a specific area of the 
north Pacific ocean. A subroutine for AP2 vessels is provided 
which can be substituted for the AP3 subroutine in the program. 
The program can be adapted to routes in other ocean areas of the 
northern hemisphere by changing the Fortran statements on a very 
small number of cards listed in the Appendix. Some suggestions 
are made for improving the program. 

II . Use of Long-Range Weather Forecasts in Ship-Routing . 

1 . The scheme of incorporating weather forecasts into the pre- 
paration of the minimal-time ship routes of this report consists 
of the following parts: 

a) The Fleet Numerical Weather Facility prepares operational 
wave predictions for periods up to 48 hours. These predict- 
ions for 1 2Z and 24Z of 26 and 27 July 1966 were used for 
the first two days in the ship voyage example of this report. 

b) The Veather Bureau's 5-day surface pressure forecast was used 
for the next three days. This forecast, which is prepared 
every Monday, Wednesday and Friday, consists of one sea-level 
pressure map per day at 1230Z. Lieutenant D. M. Truax USN [ 3 ] 
used these maps, in a Master of Science thesis project super- 
vised by Prof. Haltiner, to determine the surface winds and, 
in turn, to calculate the height, period and direction of the 
wind waves and swell. This calculated data was used at 1 2Z on 
July 28, 29 and 30 of the ship voyage example of this report. 

c) The Weather Bureau 30-day forecast was utilized for the re- 
mainder of the ship's voyage. Although not published for use 
outside the Weather Bureau, a copy of the 30-day predicted 
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moan sea-level pressure map, centered at the middle of the 
month, was mailed to Prof. Haltiner for the ship routing ex- 
periment. Lieutenant Truax [3] estimated surface vinds from 
this single map, and again the wave conditions were calcula- 
ted. These calculations were repeated on a daily basis using 
the same map for the remainder of the total forecast period. 
Since the same winds are used repeatedly during the latter 
part of the period, the forecast waves reach a steady state 
within a few days. This steady state forecast was used at 
1 2Z of July 31 and all subsequent days of the ship voyage 
example of this report. 

2. The predicted 30-day mean pressure chart has relatively weak 
pressure gradients, as would be expected from the averaging pro- 
cess. In contrast, the individual daily charts which make up 
such a mean have strong gradients in general, particularly in 
the vicinity of the migratory cyclones or low pressure areas. 
These systems have strong winds and high seas associated with 
them, which are reflected in the 30-day mean only in a very lim- 
ited fashion. The forecast procedure outlined in paragraph 1 did, 
however, show considerable skill over the use of long term month- 
ly mean charts. Nevertheless it is desirable to seek additional 
ways, possibly more accurate, of providing wave estimates for the 
latter part of a voyage extending beyond a 5-day period. One pos- 
sibility, which appears to have promise, is to develop a wave 
climatology. This could consist of utilizing the wave analyses ' 
now being prepared daily at FNVF to compute mean wave height, di- 
rection and period as a function of latitude and longitude for 
each month of the year. These data could then be compared with 
those derived from the Weather Bureau 30-day sea level pressure 
forecasts in order to ascertain the best, source of wave data for 
trans-oceanic ship routing. Such a wave climatology would have 
other applications in Naval operations as well. A further refine- 
ment in the development of a suitable wave climatology for use in 
ship routing might consist of the preparation of mean wave char- 
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acteristics not only as a function of latitude, longitude and 
month, but also separated according to weather type. The latter 
are determined largely according to the main storm tracks which 
vary from week to week as well as with seasons. Such a climato- 
logy would obviously take more effort to prepare, but would be 
a very valuable aid in ship routing. 

3. Finally, it should be mentioned that a number of groups are 
experimenting with long-range weather prediction by numerical 
integration of the hydrodynamical equations. It is expected 
that eventually such predictions will show skill for perhaps 
several weeks, and thus day-by-day wave forecasts could be made 
available for the entire period of a trans-oceanic voyage. 

III. Input and Output of the Computer Program . 

1. The program VC2AP3 , page 20 of the Appendix, was written for 
the Control Data Corporation 1604 computer in Fortran 1963, which 
is their version of the IBM Fortran IV. The 63x63 grid of the 
northern hemisphere stereographic projection of the Fleet Numer- 
ical Weather Facility was used to specify the location of the 
ocean wave data as explained in [l ] and [2]. It was desired that 
the 32764-word memory of the CDC 1604 contain wave data for as 
large a part of the 63x63 grid as possible in order to take ad- 
vantage of high-speed random access to core memory. This was ac- 
complished by packing the floating-point wave height H and wave 
direction K at a grid point into a single computer word, with 
height in the upper half and direction in the lower half word. 
Program TAPE, page 33 of the Appendix, describes how the fixed- 
point FNWF wave field data was transformed to packed floating- 
point form. The conversion of the right-shifted FNWF fixed-point 
data to floating-point was accomplished by hardware features pe- 
culiar to CDC computers, and probably not existing in IBM com- 
puters, involving normalization by addition to fixed-point octal 
2000000000000000 followed by addition to floating-point zero. 

The packing and subsequent unpacking was accomplished by CDC 
Fortran 1963 operations. The right or left shifts involved were 
simulated by division or multiplication by the integer 16777216 
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decimal. The packing and unpacking also used a CDC Fortran 1963 
masking operation which does not exist in IBM Fortran IV. 

2. The grid wave data magnetic tape output of program TAPE is 
read as input from Logical Tape Unit 1 by the main program 
VC2AP3 . It is stored in core memory as the three-dimensional 
MHD array of DIMENSION (18,32,8) with a total of 4608 words. 

The dimensions 18 and 32 correspond to FNWF stereographic pro- 
jection plane grid point indices of 8<i<25 in the direction of 
the 1 0E longitude meridian and 1 6<j<47 in the direction of the 
1 00E longitude meridian. The dimension of 8 corresponds to the 
time series of predicted wave fields as described in Section 
II, page 3 , for 1 2Z and 24Z of 26 and 27 July, 1 2Z on July 28, 

29 and 30, and the final steady-state forecast used at 1 2Z of 
July 31 and all subsequent days. The 18 by 32 rectangular field 
of data is shown on the map of Figure 1 , page 8 . The ship rou- 
ting program will not work unless all points of the ship route, 
including the initial and terminal points, are within a smaller 
16 by 30 rectangle also shown in Figure 1 , defined by 9<i<24 
and 17<j<46. A local coordinate system for the rectangular ar- 
ray is set up in the VC2AP3 program with the origin 0 at i=7 
and j=15, with the Ox and Oy axes in the direction of increas- 
ing i and j respectively. The smallest values of x and y, cor- 
responding to i=8 and j=16, are therefore x=1 and y=1 . 

3 . The MHD array of wave data in core memory may be extended 
from the present 4608 words to a maximum of 21 1 56 words to ac- 
commodate a larger ocean area or a longer time series, or both. 
This is accomplished by the relocation of COMMON feature of 
Fortran 1963 by using the following cards after the MCS con- 
trol card: 

-BINARY, 56. (The dash in column 1 is a 7 and 9 ) 

-REL0C0M. (The dash in column 1 is a zero, 7, 9 and minus ) 

-FTN,L,A,E. (The dash in column 1 is a 7 and 9 ) 

and replacing the -EXECUTE, card by -EXECUTER. A change in the 

present geometric dimensions and/or origin of the MHD array must 
be accompanied by rather obvious changes in certain cards of pro- 
grams TAPE and VC2AP3 . A list of the cards requiring a change is 
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given on page 36 of the Appendix. A change in the time dimension 
of the MHD array, corresponding to an increase in the length of 
the time series, must be accompanied by rather obvious changes 
in cards 382, 394 and 399 of subroutine TERP , with card 404 be- 
ing changed to IF ( L— 5 ) 14,5,5. An increase in the time dimen- 
sion requires changes in cards 13 and 85 of program TAPE, and in 
the format statements of cards 18, 21 , 52, 55, 64 and 67. 

4. The first two of the BCD punched card input to the program 
VC2AP3 , following the -EXECUTE, card, contain the data: 

Card 1: First line of the TI=IT title in format (6A8/) for the 

map produced by subroutine DRAW in Statement 1 1 . See ex- 
ample on page 32 of the Appendix, and explanation of the 
DRAW subroutine on page 37 of the Appendix. 

Card 2: Format (8A8,I3) of which 6A8 is the second line of the 

map title TI . The remaining part of the format is A8 for 
the DATE=KATE of the routing computation, A8 for eight 
blank Hollerith characters for a null label AL=LA on 
the map grid plot, and 13 for the NST total number of 
ships to be routed by the program. The DATE of the rout- 
ing computation corresponds to the 1 2Z hour of the first 
member of the time series in the MHD array. See example 
on page 32 of the Appendix. 

Following these two BCD cards there are groups of either 6 cards 
or one card for each ship routed by the program, depending on 
whether or not the option to plot an earlier route of a partic- 
ular ship is elected. 

Card 3: Format ( A4 , 2A8 ,F3 . 0 ,F6 . 1 ,F5.1 ,F6.1 ,F5.1 ,F3.0,I1 ,212,11 , 
F8.5,F6.3) with example on page 32 of the Appendix. The 
first A4 is the GL=LG ship identification number with 
column 1 of the card blank, used by the DRAW subroutine 
to label the terminal point of a ship route. The first 
A8 is the DATEXhKATEX date on which the ship leaves the 
initial point of its route. The second A8 is the FL=LF 
Julian date of departure, with blanks in columns 17 to 
20 inclusive, used by the DRAW subroutine to label the 
initial point of a ship route. The F3.0 is the HR hour 
of ship departure from its initial point measured from 
1 2Z on the DATE of routing, i.e. from the 1 2Z hour of 
the first member of the time series in the MHD array. 

The F6.1, F5.1, F6.1, F5.1 formats are the longitudes 
and latitudes of the initial and terminal points of the 
route, XLG1 , XLT1 , XLG2 , XLT2, with the longitudes con- 
sidered positive if east of the Greenwich meridian. The 
F3.0 format provides for the RMUL convergence factor 
discussed later on pages 11 and 16. Note the absence of 
a decimal point in the example on page 32 . The II for- 
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FIG. 1. (Left side). J207 day route of ship 666. 
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FIG. 1. (Right side). J207 day route of ship 666. 

-9- 



mat provides for NN which is either 1 or zero according 
as the option to plot an earlier route of the ship is or 
is not elected. The first of the 212 format is for the 
NSTEP reciprocal of the time step used in the integra- 
tion process, with 24 of the example indicating a time 
step of l/24 of a day. The second 12 format specifies 
the number LMAX of iterations allowed in determining the 
minimal-time route. The remaining formats II , F8.5, F6.3 
provide for the variables NP , PALF and PT used in an 
option described later in Section IV, page 13. 

Card 3 is followed by five BCD cards punched out by the state- 
ments on cards 334 to 338 of an earlier use of program VC2AP3 if 
the option NN=1 to plot an earlier track of the ship has been 
elected. If NN=0 on card 3, the remaining BCD cards of the input 
deck refer to other ships to be routed. 

5. The output of the program is a map grid for each vessel rou- 
ted, shown in Figure 1 , produced by the CALL DRAW of statement 
11. The details of subroutine DRAW are given on page 37 of the 
Appendix. If the option NN=1 has been elected, statements 1 6 to 
18 cause an earlier route of the ship to be plotted as in Figure 
2, using plus signs for daily positions and an identifying Jul- 
ian day mark for the initial point. Statements 19 to 44 cause a 
geodesic route for the ship to be computed and plotted as a so- 
lid line. Figures 1 and 2 give examples of this with the ship's 
angle of departure ALF , measured counterclockwise from the Ox 
axis, also indicated. The terminal point of the geodesic route 
plot is marked by the GL=LG identification number of the ship. 
The purpose of the geodesic route computation is to find first 
approximations to the time T and angle of departure ALF used in 
the LMAX iterations toward a minimal-time route of statements 
45 to 69. The computation of the route is abandoned if any point 
of the geodesic route falls outside of the rectangle 9<i<24 and 
17<j<46, but the geodesic route within this rectangle is plotted 
on the map grid. The minimal-time route computation is initiated 
by statement 45 only if the entire geodesic route has been com- 
puted successfully. The format of statement 71 is printed if the 
LMAX iterations result in a ship route terminal point more than 
100 nautical miles from the desired destination, together with 
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advice about how to improve the convergence of the iterative 
process. Experience to date on trans-Pacific routes indicates 
that it is desireable to use LMAX=7 and a convergence factor, 
discussed in this paragraph and later in Section IV, RMUL=10. 

The statement of card 313 gives the DIFA change in ALF for the 

I 

next of the LMAX iterations as computed from the Newton-Raphson 
equations (6) on page 19. This value of DIFA may be grossly in 
error if non-linearities in the ocean wave field happen to be 
more important than the linear terms on which the Newton-Raphson 1 
equations are based. The statement of card 315 attempts to con- 
trol the bad effects of such large non-linearities by dividing 
DIFA by a convergence factor RMUL before accepting it for the 
next iteration. If RMUL is chosen too large the convergence may 
be rapid, but the search made for a minimal-time route may be 
over such a small area of the ocean that the minimum obtained is 
local rather than global. The choice of too large a value of 
RMUL may actually lead in some cases to a local extremal with 
time greater than the geodesic route time. On the other hand 
too small a value for RMUL may lead to divergent oscillations. 
Another convergence difficulty, associated with the value of ALF 
generated by the geodesic track computation, is discussed in 
Section IY, page 13. Revisions of the program to overcome these 
convergence difficulties are suggested in Section V, page 16. 

Some monitoring of the program output will be required until the 
suggested program revisions have been accomplished. 

6. The tabulated daily ship position, wave height and direction 
for the last of the LMAX iterations are printed under the format r 
of statements 73 and 74. An example of this output is given on 
page 12, corresponding to the input IBM card on the bottom of 
page 32. Corresponding to this example the following five IBM 

cards are punched 
666 J207 13 

15.4014.3513.4312.6511.9111.4511.0810.6810.6710.35 666 J207 1 

10.3010.3210.53 0 00 0 0 0 0 666 J207 

3.05 5.06 7.09 9 . 1 6 1 0 . 9 3 1 2 . 59 1 4 . 301 6 . 141 8 . 00 1 9 . 93 J207 666 

21.8724.1326.14 000000 0J207 666 

under the formats of statements 76 and 78, which may be used for 
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PRINT output for J207 day route of ship 666 . 

TCTAL NUMBER OF VC2AP3 SHIPS ROUTED = 1 

CN JUL26 » 66 FROM 12Z TO 24Z, AND CN FOLLOWING DAY FROM 00Z TO 11Z 

ROUTE OF SHIP 666 BEGINS ON JUL26,66, JULIAN DATE = J207, 

0 HOURS AFTER I200Z CN JUL26.66 _ 

FROM LONGITUDE = -122.5 AND LATITUDE = 37.9 

TO LONGITUDE = 139.6 AND LATITUDE = 35.6 



R M U L = 10 


LMAX= 7 


NSTEP=24 


NN=0 __ 


NP = 0 






L N 1 


ALF 


T 


X { N1 ) 


XFIN 


Y ( N 1 ) 


YFIN 


0 288 


1.81954 


11.949 


13.778 


13.778 


28.357 


28.357 


1 288 


1.81954 


11.949 


14.613 


13.778 


28.159 


28.357 




1 .’82360 


11.944 










2 288 


1.82360 


11.544 


14.254 


13.778 


2 8 . 2 1 2 


28.357 




1.82600 


11.953 










3 288 


1.82600 


11.953 


14.054 


13.778 


28.311 


28.357 




1.82731 


11.930 










4 288 


1.82731 


.. .. 11.93C 


13.782 


—13.778 _J 


28.368 


28.357 




1.82731 


11.924 










5 288 


1.82731 


11.924 


13.778 


13.778 


28.357 


28.357 




1.82731 


11.924 










6 288 


1.82731 


11.924 


13.77e 


13.778 


28.357 


28.357 




1.82731 


11.924 










7 288 


1.82731 


11.924 


13.778 


13.778 


28.357 


28.357 




1.82731 


11.924 










DAYS 


LONG I 


LAT I- 


WAVE 


WAVE DIRECTION 




OF TRAVEL 


-TUDE 


TUDE 


HEIGHT 


FROM NORTH 




0 


-122.5 


'37.9 


6 


354 






1.00 


-130.3 


41.1 


5 


1 1 






2.00 


-138.9 


43.8 


2 


191 






3.00 


-148.2 


46.0 


5 


156 






4.00 


-156.9 


46.7 


21 


177 






5.00 


-164.5 


47.3 


21 


165 






6.00 


-173.1 


47.3 


21 


164 






7.00 


178.2 _ 


46.5 


18 . .. 


122 




* 


8.00 


169.7 


46.1 


19 


146 






9.00 


161.6 


43.8 


13 


124 






10.00 


154. C 


41.6 


17 


’ 357 






11.00 


146. 1_. 


38.5 _ 


10 


116 






11.92 


139.6 


35.6 


-1 


2 







GRAPH TITLED 

JOB 0574 BLEICK 

VC2AP3 DECEMBER 

HAS BEEN PLOTTED. 



BOX 6 
1966 
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some later plot of the track if the NN=1 option of card 33 of 
the program is elected in a later routing of the ship. The card 
339 statement causes the daily track positions to be plotted on 
the map grid by diamond-shaped symbols, as in Figure 1 , with the 
LFsFL Julian day identification of the initial point . Statement 
80 continues the M=1 ,NST loop to proceed with the routing of the 
next ship. 

IV. Example of Options NN=1 and NP— 1 . 

An example of the use of the NN=1 option is given in Figure 2 
where the Julian J207 day route of ship 666 is plotted with plus 
signs. The ship had departed considerably from its J207 routing 
during the first 24 hours and required a new route from longi- 
tude -130.0 and latitude 35.0 starting at 1 2Z on Julian day J208. 
It was necessary to resort to the simulation of putting HR=24 on 
this new route since a new time series MHD array starting at 1 2Z 
on J208 was not available. Considerable difficulties were encoun- 
tered in computing a global-extremum J208 day minimal-time route. 
The angle of departure ALF generated by the geodesic route com- 
putation may have bias to the extent that successive track iter- 
ations always go around a storm area on the same side, when actu- 
ally the other side would be a better choice. The geodesic track 
time was T=1 1.385 days and angle of departure ALF=1. 73554 radians 
This ALF was a poor initial approximation to that required for th 
final successful minimal-time route of Run 3 , as indicated by the 
successive runs of the following table: 



Type Graph Converged 

of mode NP=1 input VC2AP3 output • 



Run 


NP 


Route 


Fig. 2 


RMUL 


LMAX 


PALF 


FT 


ALF 


T 


0 

1 


0 


Geo- 

desic 

Local 

extreme 


Solid 

line 

Dashed 

line 


10 


1 0 






1 .73554 
1 .749 


1 1 .385 
11.567 


2 


1 


Local 

extreme 


Dashed 

line 


4 


6 


1 .81 200 


10.967 


1 .814 


1 1 .682 


3 


1 


Global 

extreme 


Diamond 

symbols 


999 


7 


1 .90000 


1 1 .250 


1 .90003 


11.264 
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FIG. 2. (Right side). J208 day route with options NN=1 & NP=1 
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The program VC2AP3 provides the option to override the values of 
ALF and T generated by the geodesic route computation by using 
NP=1 in column 54 of the first ship-related IBM input card des- 
cribed on page 7, and to start the LMAX iterations with the val- 
ues ALF=PALF (format F8.5) and T=PT (format F6.3) from columns 55 
to 68 inclusive of this IBM card. The values of PALF and PT in 
Runs 2 and 3 of the above table were chosen by inspection of the 
iterations produced in earlier runs. Note that Runs 1 and 2 con- 
verged to local time-extremal routes requiring more time than the 
geodesic route. The successful Run 3 required the input RMUL=999 
and PALF=1 .90000 radians, the former probably being somewhat lar- 
ger than really required for convergence. The nearby input RMUL= 
100 and PALF=1 .90000 were tried without obtaining convergence. 

The extreme value of the convergence factor RMUL=999 used in Run 
3 and the resulting sensitivity in the required value of the in- 
put angle of departure PALF suggest that the non-linear terms ne- 
glected in the Newton-Raphson equations are very strong on the 
J208 day route. 

V. Conclusions . 

1 . It is suggested that any future revision of program VC2AP3 
take steps to reduce the computer running time of a trans-Pacific 
track iteration from the present l! minute, and 47 seconds to an 
attainable 13 seconds. This is accomplished easily by eliminating 
the calculation of 128 sine and cosine functions on each entry of 
subroutine TERP. Under this proposed scheme the MHD array would 
store a grid-point wave height H only in each word. Two addition- 
al 4608-word arrays would store the the CK and SK trigonometric 
functions, computed once only in a modified program TAPE. Assum- 
ing that the time dimension of these three arrays remains at 8, 
the CDC 1604 memory capacity would even allow the geometric di- 
mensions to be chosen to cover 881 grid points, an increase of 

/ 

about 53$ from the present 576 grid points of the MHD array. This 
proposal would eliminate the normalizing and masking hardware fea- 
tures of the program peculiar to CDC computers, and make the pro- 
gram useable in other Fortran IV systems. Even larger dimensions 



for the three arrays could be achieved by using a more realistic 
graph plotting program to replace subroutine DRAW with its exhor— 
bitant demand for 5760 words of core memory. Care would have to 
be exercised in implementing this proposal since the X(900) and 
Y(900) arrays of subroutine DRAW are used by program VC2AP3 for 
other memory-conserving purposes. 

2. It is suggested that any future revision of program VC2AP3 » 
try to go beyond the present linear Newton-Raphson scheme, and 
attempt to include quadratic terms also. Such a change might 
overcome the convergence difficulties disclosed in Section IV. 

If such a quadratic scheme were found and used successfully, it 
may be possible to obtain a considerable reduction in the LMAX 
number of track iterations required to achieve convergence. 

3. In order to reduce the amount of monitoring of output re- 
quired in program VC2AP3 , it is suggested that any future revi- 
sion of the program attempt to include as far as possible the 

following recommendations of Schmieg, [5] pages 30, 31 and 37: 

. conditions 

a) Check the Legendre and Weierstras s A f or a minimal-time route 

at each time step of the numerical integration. This check 
may give an automatic method of eliminating local time- 
extremal routes which are not the desired global minimal- 
time route. 

b) Examine the magnitude of the Hamiltonian at the end of each 
of the LMAX iterations to determine whether it has diminish- 

* 

ed from its value at the end of the previous iteration. If 
not diminished choose an appropriate convergence factor 
RMUL to get ALF for the next track iteration. 
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VII. Mathematical Theory . 

1. A misprint in equation ( 1 6) of [l] and [2] requires correc- 

2 

tion by replacing R by R . The 6p of this equation is an approx- 
imation in that certain terms obtained in the differentiation of 
equation ( 1 0 ) have been dropped as being unimportant. Vhile the 
approximation is useful for a short ship route, it is desireable 
to use the following completely accurate expression for a long 
trans-Pacific route 



6p = ^[R 2 |E|A" 2 6a + (V V-V V )6x + (V V-V V )6y], 

r D L I*-.!" PX p X DV T> V J J 



py 



p y 



(D 



2 2 

where D=V +2V -W . Use of this new expression for 6p with 
P PP r * 

equations (14), (16) and (18) of [l] and [2] gives the following 

pair of inhomogeneous differential equations as a replacement 

for equation ( 19 ) of these references 



2 j 

^(X,6x +Ul 6 y )=^[x(VY py -V p V y )-p(W px -V i> V x )](X,6xH- tll 6y)- s inaJf- (|) 6a 

( 2 ) 

2 3 

^(X 2 6x +l r 2 6y)=^[x(VY -V V )-M(W -V V x )](x 2 6x +t , 2 6y) + cosaJfL (5) 6a 



Let yj be a solution of either of (2), rendered homogeneous by 
putting 6a=0 , and with y^(0)=1 . Then a solution of (2) with 
6a^0, and with vanishing initial values is 



where 



(\ 1 6 x+ 4 1 6 y) T = -y 7 (T)y g (T) sin a 6 a 
(\ 2 (Sx +U 2 ( 5 y)x = +y 7 (T)y g (T) cos a 6 a 



y 8 (T) 



J °T 



^(R/A) J dt. 



(3) 

(4) 
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The solution of (3) gives the following new version of equa- 
tions (19) of [l] and [2] 

[6x,6y] T = y 7 (T)y 8 (T)6a[-u,X] T /|E(T)| , ( 5 ) 



and the Newton-Raphson equations (22) of these references 
replaced by 

x ( T ) AT - [My :7 y 8 /|E| ] T 6a = Ax(T) 
y(T) AT + [Xy 7 y g /|E| ] T 6a = Ay(T) . 



are 

( 6 ) 



2. The subroutines AP3 and AP2 of the Appendix are based on 
the work of James [ 4 ]. His speed versus wave height data for 
these vessels have been approximated by arcs of hyperbolas as 
explained in [l] and [ 2 ]. 
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IX. Appendix . 

—COOP » BOX 6, BLEICK*I/l/0/49/S/56/57/E/45»54*21 *10000.0. VC2AP3 - 6 DEC 66. 
-FTN »L »A»E. 

PROGRAM VC2AP3 0 



C YVARS ( I ) » LAMBDA 1 YVARS(2)=MU1 YVARS < 3 ) *LAMBDA2 YVARS(4)=MU2 
C YVARS ( 5 ) *X YVARS ( 6 ) =Y YVARS(7)«S OR Y7 YVARS(8)»Y8 

DIMENSION MHD( 18 » 32 * 8 ) * X ( 900 ) » Y( 900 ) *RX ( 10*90) *RY ( 10 *90 ) * 1 

+ IT(12)*TI(12)*C(4) * YVARS ( 8 ) » YC ( 8 ) *AK (4*8) * DY ( 8 ) 2 

COMMON MHD»TC»YC*H*HX*HY.CK,CKX.CKY.SK*SKX*SKY*XK*YVARS*XLG»XLT* 3 
+ A, B » CC * DA *DB » DC *LR 4 

EQUIVALENCE ( IT,TI)»(LA»AL)»( (CATE* DATE) • { LP.PL) * ( LG.GL) * ( LF.FL) • 5 

+ (X.RX) *(Y*RY) * (KATEX.DATEX) 6 

REWIND 1 7 

C(l) = 0.0 8 

C ( 2 ) = 0.5 9 

C ( 3 ) = 0.5 10 

C ( 4 ) = 1.0 11 

C READ MAP GRID DATA FOR DRAW SUBROUTINE* AND WAVE FIELD MATRIX 

READ ( 1 ) X » Y 12 

READ < 1 ) MHD 12. 

C READ MAP TITLE* DATE OF ROUTING COMPUTATION* MAP GRID PLOT LABEL* 

C AND TOTAL NUMBER OF SHIPS ROUTED 

READ (50*1 ) T I *DATE*AL*NST 13 

1 FORMAT (6 A8/8A8 *13) 14 

WR I TE ( 51 * 2 ) NST *KATE 15 

2 FORMAT ( 39H1TOTAL NUMBER OF VC2AP3 SHIPS ROUTED “ I3/1X3HON A8*54H 16 

+FROM 12Z TO 24Z* AND ON FOLLOWING DAY <FROM OOZ TO 11Z/) 17 

REWIND 1 18 

DO 80 M»1 * NST 19 

IF ( M-l ) 10*11*10 20 

C READ MAP GRID DATA FOR DRAW SUBROUTINE 

10 READ ( 1 ) X* Y 21 

REWIND 1 22 

C DRAW MAP GRID 

11 CALL DRAW ( 386 ,X * Y * 1 *0 *L A * I T » 2 • ♦ 2 • *0 *0 * 2 * 2 *9 » 1 5 ,0 , LAST ) 23 



C READ SHIP IDENTIFICATION NUMBER* DATE AND HOUR OF DEPARTURE* COORDINATES 
C OF TRACK. END POINTS. CONVERGENCE FACTOR* OPTION TO PLOT EARLIER 

C TRACK, TIME STEP RECIPROCAL, AND NUMBER OF ITERATIONS 

READ <50, 14) GL.DATEX*FL*HR*XLG1*XLT1*XLG2*XLT2*RMUL.NN,NSTEP*LMAX*24 
+ NP,PALF,PT 24, 

14 FORMAT <A4, 2A8 * F3 .0 * F6. 1 *F5. 1 , F6. 1 *F5. 1 ♦ F3. 0*11*212. II. F8.5.F6. 3) 25 

RSTEP * NST EP 26 

WRITE! 51,15) LG ,KATEX »LF*HR*KATE»XLG1 *XLT1 *XLG2 *XLT 2 *RMUL*LMAX* 27 

+ NSTEP ,NN ,NP 28 

15 FORMAT ( 15HOROUTE OF SHIP A4*llH BEGINS ON A8*16H* JULIAN DATE * A429 

+,1H,/1XF3.0*22H HOURS AFTER 1200Z ON A8/19H FROM LONGITUDE « F6.130 
+»16H AND LATITUDE = F6.1/19H TO LONGITUDE * F6.1.16H AND LATITU31 
+DE = F6.1//6H RMUL=F5.0*3X5HLMAX*I2*3X6HNSTEP«I2*3X3HNN«I1*3X3HNP»32 
+ 11//) 32, 






- 20 - 



C CHECK ON OPTION TO PLOT EARLIER TRACK 

IF (NN) 16 ♦ 19 * 16 33 

16 READ ( 50 » 1 7 ) GL»PL»NK 34 

17 FORMAT ( 2 A8 ♦ I 2 ) 35 

READ (50*29) ( X ( I ) * 1 3 1 , 20 ) ♦ ( Y ( I ) ♦ I =1 ,20 ) 36 

29 FORMAT (10F5.2) 37 

CALL DRAW ( NK ♦ X » Y * 2 * 2 » LP » I T , 2 . , 2 . , 0 , 0 ♦ 2 , 2 , 9 , 1 5 » 0 , LAST ) 38 

WR I TE ( 51 * 18 ) LG * LP 39 

18 FORMAT! 23 HO EARLIER ROUTE OF SHIP A8,15H ON JULIAN DAY A4/66H HAS 40 
+BEEN PLOTTED USING PLUS SIGNS FOR SUCCESSIVE DAILY POSITIONS/) 41 

C COMPUTATION OF GEODESIC TRACK 

19 ARG = (XLG1-10. )/57. 29577951 '* 42 

C0SLG1 3 COSF(ARG) 43 

S I NLG1 = SINF(ARG) 44 

ARG 3 (XLG2-10. )/57. 29577951 45 

C0SLG2 3 COSF(ARG) 46 

SINLG2* S INF ( ARG ) 47 

ARG = XLT1/57. 29577951 48 

C0SLT1= COSF(ARG) 49 

S I NLT 1 3 S INF ( ARG) 50 

ARG 3 XLT2/57. 29577951 51 

COSLT2 = COSF(ARG) 52 

S I NLT2 3 SINF(ARG) 53 

EL 3 S I NLT 2*COSLT 1*S I NLG1 - COSLT2*S I NLT 1*S I NlG2 54 

EM =-SlNLT2*COSLTl*COSLGl + COSLT2*S INLT1*C0SLG2 55 

EN =(SINLG2*COSLGl-COSLG2*SINLGl )*COSLTl*COSLT2 56 

ROOT = SQRTF ( EL*EL + EM*EM + EN*EN ) 57 

EL = EL/ROOT 58 

EM = EM/ROOT . 59 

EN = EN/ROOT , 60 

PR1= 31.205*COSLTl/< 1.+SINLT1) 61 

XI = PR 1*C0SLG1 62 

Y1 = PR 1*S I NLG 1 63 

PR2= 31.205*COSLT2/( 1.+SINLT2) 64 

X2 = PR2*COSLG2 .. 65 

Y2 = PR2*SI NLG2 66 

DELX * X2 - XI ■' 67 

DELY = Y2 - Y1 68 

S12 = SQRTF ( DELX* DELX + DELY*DELY) 69 

ARC= S 12 * 70 

IF (XLG2-XLG1) 20,21*20 71 

20 ARG= ABSF ( EN/62,41 ) 72 

ARC= ASINF(ARG*S12 )/ARG 73 

21 COSA = -EN*Y1/31.205 + EM 74 

SINA = EN*X1/31.205 - EL 75 

IF (.COSA) 23 » 22 » 23 76 

22 ALF = SIGNF(1«5707963268*SINA) 7 77 

GO TO 27 78 

23 ALF = ATANF(SINA/COSA) 79 

IF (COSA) 24,27,27 80 
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l o i i 






2 A 


IF (SINA) 


26,25,25 


81 


25 


ALF = ALF + 


3.1 A 15926536 


82 




GO TO 27 




83 


26 


ALF = ALF - 


3.1415926536 


84 


27 


N3 = 0 




85 




X ( 1 ) = XI + 


24 o 


86 




Y ( 1 ) = Y 1 + 


16. 


87 




XFIN = X2 + 


24. 


88 




YFIN = Y2 + 


16. 


89 




LR = 0 




90 




STEP = 1 • O/RSTEP 


91 




TAU = 0.0 




92 




S =0.0 




93 



TVAR = HR/24. 
YVARS ( 5 ) = X ( 1 ) 
YVARS ( 6 ) = Y ( 1 ) 



0.0 



YVARS ( 7 ) 

N1 = 1 
N2 = 1 

DO AO K=2 » 900 
DO 32 1=1, A 

TC = C C I >*STEP + TVAR 
DO 31 J= 5 » 7 

31 YC(J> = C( I )*AK( 1-1, J) + YVARS(J) 

IF ( ABSF ( YC ( 5 ) - 9.5)- 7.5) 97,38,38 

97 IF (ABSF(YC(6)-16*5)-1A,5) 98,38,38 

98 CALL TERP 
CALL AP3 



DELX 

DELY 

COSP 

SINP 

COST 



YC ( 5 ) - 2A • 

YC ( 6 ) - 16. 

-DELY*EN/31.205 + EM 
DELX*EN/31 ,205 - EL 
COSP-*CK + SINP*SK 
A2MC2= A*A - CC*CC 
SPMK = (SINP*CK-COSP*SK)/B 
SI NTB= SPMK*SPMK*A2MC2 

V = A2MC2 / ( CC*COST + SQRTF ( COST *COST +SINTB ) *A ) 

EMF I = ( DELX*DELX + DELY*DELY + 973.75 ) /10A3 .6387A3 

CAPV = V* EM FI/0.5660416667 

DY ( 5 ) = CAPV*COSP 

DY ( 6 ) = CAPV #S I NP 

DY ( 7 ) = CAPV 

DO 32 J= 5 , 7 

32 A K ( I , J ) = STEP*DY( J) 

DO 33 J= 5 » 7 

33 YVARS ( J ) = ( AK(1,J)+2.*AK(2>J)+2.*AK<3»J)+AK(4,J) )/6, 

TVAR = TVAR + STEP / 

X ( K ) = YVARS ( 5 ) 

Y(K) = YVARS ( 6 ) 

Ml = K 



9 A 

95 

96 

97 

98 

99 

100 
101 
102 
103 
10A 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 
117 

lie 

ns 

120 

121 

122 

12 : 

12 L 

12? 

YVARS ( J ) 126 

121 

121 
1 2 
1 3 ( 
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34 



N2 = K 131 

IF ( YVARS ( 7 ) “ARC ) 35,34,34 132 

RAT = { ARC-S ) / ( YVARS ( 7 ) -S ) 133 

T = STEP*RAT + TAU 133.1 

X ( K ) = <X!K)-X(K-1 ) )*RAT + X(K-l) 133.2 

Y ( K ) = (Y( K)-Y( K-l ) )*RAT + Y(K-l) 133.3 

N2 = K+l 134 

X(N2) = XFIN 135 

Y(N2) = YFIN , 136 

GO TO 41 137 

35 S = YVARS ( 7 ) * 138 

TAU = TAU + STEP 139 

IF ( K-900 ) 36,38,38 140 

36 IF ( ABSF ( X ( K ) - 9.5)- 7.5) 37,38,38 141 

37 IF ( ABSF(Y( KJ-16.5 )-14.5) 40,38,38 142 

38 T = TAU 143 

WR I TE ( 5 1 , 39 ) LG 144 

39 FORMAT (61HOMORE THAN 899 INTEGRATION STEPS OR WAVE DATA FIELD EXCE145 

+EDED./21H OTS ROUTING OF SHIP A4»4X39HABANDONED BUT GEODESIC TRACK146 
+ IS PLOTTED/) 147 

N3 = 1 148 

GO TO 41 ' 149 

40 CONTINUE 150 

41 L = 0 151 

WR I TE ( 5 1 , 42 ) 152 

42 FORMAT! 4X 1HL4X2HN16X3HALF7X1HT7X5HX < Nl) 4X4HX F I N 5X5HY ( N1 J4X4HYF IN/ ) 153 
PRINT WEIGHTING FACTOR ALPHA AND TIME T OF GEODESIC TRACK 

WR I TE ( 5 1 ,43 ) L, Nl * ALF » T ,X < N 1 ) , XFIN, Y(N1) , YFIN 154 

43 FORMAT ( I 5 , 1 6 , FI 1 . 5 , 5F9 . 3 ) 155 

ROTATE AND TRANSLATE AXES TO PLOT GEODESIC TRACK ON MAP GRID 

DO 44 1=1, N2 156 

TEMP = .97780241408*X< I ) - • 20952908873*Y( I ) + 3.0032502718 157 

Y ( I ) = .20952908873*X( I ) + .97780241408*Y ( I ) - 4.4699538929 158 

44 X ( I ) = TEMP 159 

CALL DRAW ( N2 , X , Y , N3+2 , 0 , LG , I T , 2 . , 2 . ♦ 0 , 0 , 2 , 2 , 9 , 1 5 , 0 ♦ LAST ) 160 

IF ( N3 ) 80,45,80 161 

PREPARE FOR LMAX ITERATIONS TOWARD MINIMAL-TIME TRACK 

45 TC = HR/24. 162 

X<1) = XI + 24. 163 

Y ( 1 ) = Y1 + 16. 164 

YC ( 5 ) = X ( 1 ) 165 

YC(6> = Y ( 1 ) 166 

CALL TERP 167 

DO 9 1=2,399 •' 168 

X! I ) = 0.0 x 169 

9 Y ( I ) = 0.0 , 170 

X ( 1 0 1 ) = H / 171 

YVARS! 5 ) = X<1) 172 

YVARS ( 6 ) = Y ( 1 ) 173 

CALL ANGLE 174 
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Y < 101 ) « XK 175 

X ( 201 ) = XLG1 176 

Y ( 201 ) * XLT1 177 

LR » 1 178 

IF <NP) 81,32,81 178.1 

81 ALF 3 PALF 178.2 

T - PT 178.3 

COSA * COSF(ALF) 178.4 

SINA * S I NF ( ALF ) 178.5 

82 DO 69 L 3 1 » LMAX 179 

TVAR = HR/24. „ 180 

TAU 3 0.0 181 

N1 3 XI NTF ( RSTEP * T) 182 

XN1 3 N1 183 

STEP1 3 1.0/RSTEP 184 

FSTEP ■ -XN1/R5TEP H T 185 

N1 3 N1 + 2 186 

DO 46 1*1,8 187 

46 YVARS( I ) * 0.0 % 188 

YVARS(l) * 1.0 _ 189 

YVARS ( 4 ) 3 1.0 * 190 

YVARS ( 5 ) * XQ) 191 

YVARS ( 6 ) * Y ( 1 ) 192 

YVARS ( 7 ) 3 1.0 192.1 

NX 3 1 193 

DO 66 K=2,N1 194 

STEP 3 STEP1 195 

IF ( K— N 1 ) 48,47,40 196 

47 STEP 3 FSTEP 197 

48 DO 52 1*1,4 198 

TC * C(I)*STEP + TVAR 199 

DO 49 J*l.,8 200 

49 YC(J) 3 C( I )*AK( 1-1 » J) + YVARS(J) 201 

IF ( A8SF ( YC ( 5 ) - 9.5)- 7.5) 99,65,65 . 202 

99 IF (ABSF(YC(6)-16. 51-14.5) 100,65,65 203 

100 XLAM 3 YC(l)*COSA + YC(3)*SINA 204 

XMU = YC ( 2 ) *COSA + YC(4)*SINA 205 

CLAM 3 SQRTF ( XLAM*XLAM + XMU*XMU ) 206 

CALL TERP 207 

CALL AP3 208 

DKX 3 CK*SKX - SK*CXX 209 

DKY 3 CK*SKY - SK#CKY 210 

ABS 3 (XLAM*CK + XMU*SK ) *A/CLAM 211 

ORD 3 (XMU *CK ~XLAM*SK)*D/CLAM 212 

HYP 3 SQRTF (ABS*ABS + ORD*QRD) 213 

SINB* ORD/HYP 214 

COSB 3 ABS/ HYP ' 215 

VMAJ= A * COSB - CC 216 

VMI N 3 B * SINB 217 
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V = SQRTF ( VMAJ*VMA J + VMIN*VMIN) 210 

COSP* ( CK.*VMAJ - SK*VM IN ) /V 219 

SINP* (SK*VMAJ + CK.*VM IN ) /V 220 

COST® VMAJ/V 221 

VBR = V/8. 5660416667 222 

DELX = YC ( 5 > - 24* 223 

DELY* YC ( 6 ) - 16* 224 

EMFI * ( DELX*DELX + DELY*DELY ♦ 973.75) /1043.638743 225 

CAP V = VBR * EMFI 226 

DY < 5 ) = CAPV * COSP 227 

DY ( 6 ) = CAPV * SINP 228 

EMF I X= DELX/521. 8193715 229 

EMF I Y= DELY/521. 8193715 230 

B2MA2* B*B - A*A * 231 

AMCCB* -CC*COS6 + A 232 

ACPOCB = A*CC + B2MA2*COSB 232.1 

RAT * (ACPDCB*SINB/B) /AMCCB 233 

VP » RAT * V 234 

CAP VP* RAT * CAPV 235 

QUO = V/AMCCB 236 

DIV = 1.0 237 

IF (RAT) 12.13.1 2 238 

01 V = RAT*RAT + ( B2MA2*SINB*SINB/B - VP*COST/SINB)*QUO*QUO/B ♦ 1. 239 
DET = YC( 1 ) *YC( 4) - YC(2)*YC(3) 240 

QUO = SQRTF ( RAT*RAT + 1.0)/CLAM 241 

FORCE* CAPV*DET*DET*QUO*QUO*QUO/DIV 242 

FNUM = -RAT * AMCCB 243 

OBOH = ( DA*COSB - DC)*COSB + A*DB*SINB*S INB/B 244 

RATX = ( FNUM*DKX + HX*DBDH ) /AMCCB 245 

RATY * ( FNUM*DKY + HY*DBDH ) /AMCCB 246 

CAPVX* ( RATX*EMFI + EMFIX)*V8R 247 

CAPVY* <RATY*EMF1 + EMFI Y ) *VBR 248 

Fl=( ( ( ( <A+AMCC8)*B2MA2*C0SB+A*A*CC)*C0$B-(CC*CC+62MA2)*A) /AMCCB)/ 248.01 
+ AMCCB) /B 248.02 

F2 * ( <DA/A-DB/B)*COSB - DC/A) #SINB*F1 248.03 

F3 = (CC*SINB*F1/A)/CAPV 248.04 

F4 = ( (ACPOCB*COSB - B*B)*F1/A)/B 240.05 

F5 = ( ( (B*DB-A*DA)*2.*C0SB-k:C*DA+A*DC)/ACPDCB-D8/B-MDC*C0SB-£>A)/ 248.06 

+ AMCCB ) *RAT 248.07 

QUOX = ( F 2+F5 ) *HX + F3*CAPVX + F4*DKX 248.08 

QUOY = ( F 2+F5 ) *HY + F3»CAPVY + F4*DKY 248.09 

CAPVPX * CAPVX*RAT + CAPV*QUOX 248.10 

CAPVPY = CAPVY*RAT + CAPV*QUOY 248.11 

F6 = ( -RAT*CAPVX + CAPVPX )*QUO/DIV 248.12 

F7 * ( -RAT*CAPVY + CAPVPY »*QUO/DIV 248.13 

F8 = F7*XLAM - F6*XMU 248.14 

DY ( 7 ) = F8*YC ( 7 ) 248.15 

DY ( 8 ) * FORCE/YC( 7 ) . . . 248.16 
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SUM » YC ( 1 ) *COSP + YC (2 )*SINp 249 

DY ( 1 )* -CAP VX * SUM 250 

DY C 2 ) » -CAPVY * SUM 251 

SUM « YC ( 3 ) *COSP + YC(4)*SINP 252 

DY ( 3 ) * -CAPVX * SUM 253 

DY ( 4 )* -CAPVY * SUM 254 

DO 52 J*l,8 255 

52 AX ( I » J ) = STEP * DY < J ) 256 

DO 53 J B 1 » 8 257 

53 YVARS(J) = ( AX< 1 ♦ J ) +2 ♦ *AX < 2 • J > +2.*AX 1 3 , J )*AX ! 4, J > */6« ♦ YVARS<J)258 

® TVAR = TVAR + STEP 259 

TAU = TAU + STEP 260 

IF (Nl-X) 54,56*54 261 

54 IF ( LMAX-L ) 62*55,6 2 262 

55 IF < (X-1)/NSTEP+1-NX> 62*62*56 , 263 

56 NX = NX + 1 264 

XLAM » YVARS!l)*COSA + YVARS!3)*SINA 265 

XMU = YVARS!2)*COSA + YVARS<4)*S1NA 266 

CLAM « SQRTF!XLAM*XLAM + XMU*XMU) 267 

. YC ( 5 ) = YVARS { 5 ) 268 

YC ( 6 ) = YVARS(6 ) 269 

LR » 0 270 

CALL TERP 271 

CALL AP3 272 

LR « 1 273 

ABS » < XLAM*CX + XMU*SX ) *A/CLAM 274 

ORD * ( XMU*CX - XLAM*SX)*B/CLAM 275 

HYP » SQRTF ( ABS*ABS ORD*GRD) 276 

VMAJa A * ABS/HYP - CC • . 277 

' VMIN» B * ORD/HYP ’ 278 j 

. V » SQRTF (VMAJ*VMAJ + VMIN*VMIN> 279 

COSP= ( CX*VMAJ - SX*VM IN ) /V 280 

SINP= { SX*VMAJ CK*VMIN)/V 281 

DET = YVARS ( 1 ) * YVARS! 4 ) - YVARS! 3)*YVARSI 2 > 282 



CONA»-XMU*YVARS( 7 I *YVARS! 8 >/D£T 282* 

C0NB“ XLAM* YVARS! 7 ) *YVAR5! 8 )/DET 282* 



IF (LMAX-L) 61, 60. 61* 283 

60 X ( NX ) - YVARS ( 5 ) 284 

Y ( NX ) a YVARS! 6 ) 285 

X! NK+100 ) - H 286 

X! NX+300) - TAU 287 

CALL ANCLE ' 288 
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Y ( NK+ 1 00 ) = XK 289 . 

X ( NK+200 ) = XLG 290 . 

Y ( NK+200 ) = XLT 291 

61 IF (Nl-K) 62 *67*62 292 

62 DELX = YVARS ( 5 ) - X ( 1 ) 293 

DELY = YVARS ( 6 ) - Y ( 1 ) 294 

IF ( DELX*DELX + DELY*DELY - S12*S12) 63,65,65 295 • 

63 IF <ABSF<YVARS< 5)- 9.5)- 7.5) 64,65,65 296 

64 IF (ABSF(YVARS<6)-16*5)-14.5) 66,65,65 297 

65 N1 = K 298 

T = TAU 299 ' 

GO TO 56 300 

66 CONTINUE 301 . 

PRINT ALPHA, T, X, AND Y AT END OF EACH ITERATION 

67 WR I T E < 5 1 , 43 ) L » N1 , ALF »T , YVARS < 5 ) » XF I N , YVARS < 6 ) ♦ YF IN 302 

DELX = YVARS ( 5 ) - 24. 303 

DELY = YVARS ( 6 ) - 16. 304 

EMF I = ( DELX*DELX + DELY*DELY + 973.75) /1043 .638743 305 

DIFX = XF I N - YVARS ( 5 ) 306 

DIFY = YF I N - YVARS < 6 ) 307 

CAP V = V * EMFI /8. 5660416667 308 

XDOT = CAPV * COSP 309 

YDOT = CAPV * SINP 310 

DETER= XDOT *CONB - YDOT *CONA 311 

DIFT = <CONB*DIFX - CONA*DIFY ) /DETER 312 

DIFA = ( XDOT *DI FY - YDOT*DI FX ) /DETER 313 

T = DIFT + T 314 

ALF = DIFA/RMUL + ALF 315 

COSA = COSF(ALF) 316 

SINA = SINF(ALF) , 317 

PRINT NEW VALUES OF ALPHA AND T 

WR I TE ( 5 1 , 50 ) ALF , T 318 

50 FORMAT < 1 1XF1 1 . 5 , F9. 3 ) 318.1 

69 CONTINUE 319 

IF < DI FX*D I FX + DI FY*DI FY - EMF I *EMF I * . 2366 ) . 72,72,70 320 . 

70 WR I TE ( 51,71) LG 321 

71 FORMAT ( 20H0 OTS ROUTE OF SHIP A4,47H MORE THAN 100 MILES FROM DEST322 

+ 1 NAT ION BUT TRACK/69H IS PLOTTED. INCREASE RMUL OR LMAX , OR BOTH, 323 
+ TO IMPROVE CONVERGENCE.) 324 

TABULATE FINAL TRACK DAILY POSITION, WAVE HEIGHT AND DIRECTION 

72 WR I TE ( 5 1 , 73 ) 325 

73 FORMAT ( 1 H0/4X4HDA YS7X 5HLONG I 4X5HLA T I -5X4HWA VE5X 14HWAVE DIRECTION/326 

+2X9HOF TRAVEL4X5H-TUDE4X4HTUDE5X6HHEIGHT6X10HFROM NORTH/) 327 

WRITE(51,74) ( X ( K+300 ) ,X(K+200) , Y < K+200 ) , X ( K+l 00 ) , Y ( K+100 ) , K = 1 ,NK ) 328 • 

74 FORMAT <F9.2»F11.1,F8.1,F9.0,F14.0) 329 

ROTATE AND TRANSLATE AXES FOR PLOT OF DAILY POSITIONS 

DO 75 1=1, NK 1 z 330 

TEMP = .97780241408*X< I ) - • 20952908873*Y ( I ) + 3.0032502718 331 

Y ( I ) = .209 52908873*X < I ) + . 97 7802 4 1408* Y ( I ) - 4.4699538929 332 

75 X < I ) = TEMP 333 
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c 



c 



PUNCH 11 CARDS USEABLE FOR A LATER PLOT OF TRACK 
WRI TE( 52 »76 ) LG>LF>NK 
76 FORMAT ( 2A8 . 1 2 .61X1H ) 

WRI TEC 52 >78 ) C ( (RXC 1 > J ) > I «1 > 10 ) >LG#LF. J ) > J«1 >2) > 

♦ (( (Rr<I.J>.I-l>10) >LF>LG.J) >J»1>2) 

78 FORMAT ( 10F5 • 2 > 2A8 1 1 2 > 11X1H ) 

CALL DRAW (NK>X»Y»3.4.LF.IT.2 # >2« >0>0 >2 > 2 >9> 15 1 0> LAST ) 

WRI TEC 51 >93 ) 

93 FORMAT (1H1) 

PROCEED TO COMPUTE THE ROUTE OF NEXT SHIP 
80 CONTINUE 
STOP 
END 

SUBROUTINE TERP 

DIMENSION MHD ( 4608 ) > YC ( 8 ) >HT ( 4 >4 ) >CT ( 4 >4 ) *ST ( 4>4 ) >YVARS C 8 ) >P(4) • 

+ Q< A) >PX( 4 ) >QY( A) >HDC 4 ) >CD( 4 ) >SD( A) >HS( 4 ) .CSC 4 ) >SS( A ) >HP( A ) >CP( A) 
+ SPC 4 ) >HPX( 4 ) >HPY( 4 > > CPXC 4 ) .CPYC 4 ) >SPXC 4 ) >SP Y 1 4 1 *CC 4 ) > XHO( 4608 ) 
COMMON MHD >TC>YC>H>HX>HY >CK >CKX>CKY>SK >SKX»SKY »XK . YVARS.XLG.XLT • 



334 

335 

336 

337 

338 

339 

340 

341 



342 

343 

344 
343 
346 

.34’, 

34C 

34‘ 



+ A>B>CC >DA>DB>DC >LR 

EQUIVALENCE ( MHD.XHD) > ( ID.ARG) 
MASK* 77777777B 
DTC = 2 >*TC 
L = XI NTF (DTC) 

IF (L-3) 1 >1 >7 

1 TT * ( -INTF ( DTO+DTC) *2. - 1. 
TP1* TT + 1. 

TM 1= TT - 1. 

T2M= TP1*T Ml 
IF (U 2 >2 >3 

2 K4 = 3 

TM3= TT - 3, 

CC1)» TMl*TM3/8 • 
C(2)=-TPl*TM3/4. 

C( 3 ) = T2M/8. 

GO TO 16 

3 K4 = 4 

IF ( L-2 ) 4.4»6 

4 G = (3>*TT+2. >*TT - 9. 

F = -4.*TT + G 

CC1)= -T2M*TM1/16» 

C ( 2 ) = G*TM1/16. 

5 C C 3 ) = -F*TP1/16> 

C C 4 ) =* T2M*TP1/16. 

GO TO 15 

6 C C 1 » = -T2M*TM1/16. 

C( 2 ) * ( ( 2 >*TT+1 • )*TT-7.)*TM1/12. 
C ( 3 ) = ( C l>-TT)*TT+4.)*TPl/8.' 

C( 4 ) * T2M*TPl/48. 

GO TO 15 

7 L « XINTFCTC-2.) + 4 



351 



35 

35 

35 

35 

35 

35 

35 

36 
36 
36 
36 
36 
36 
36 
36' 
36t 
3<r 
3)1 
3T I 
3 I 
33 I 
3*1 
3 5 1 
3 I 
3 I 

3 sl 

3 I 
3 I 

I 
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IF ( L— 8 ) 9 * 8 » 8 332 

8 K4 = 1 333' 

L = 7 384 

C( 1 ) = 1. 385 

GO TO 16 386 

9 TT = (-INTF(TC)+TC)*2. - 1. 387 

TP 1= TT + 1. 388. 

TM1= TT - 1. 389 

T 2M= TP 1*TM 1 390 

G = ( 3.*TT+2. )*TT - 9, 391 

F = -4.*TT + G 392 

C( 1 ) = -T2M*TM1/16. 393- 

IF ( L-7 ) 11*10*8 394 

10 K.4 = 2 395 

C ( 2 ) = ( G*TM1 + (T2M-F)*TP1 )/l6. 396 

GO TO 15 397 

11 C ( 2 ) = G*TM1 / 16 • 398' 

IF ( L-6 ) 13*12*10 399 

12 K.4 = 3 . 400 

C ( 3 ) = ( T2M-F)*TP1/16. 401 

GO TO 15 402 

13 K4 = 4 403 

IF ( L-5 ) 14*5*12 404 

14 C(l) = -T2M*TMl/6. 405 

C ( 2 ) = ( ( 5.*TT+2. )*TT-11. )*TM1/16. 406 

C ( 3 ) * ( (-5.*TT+4.)*TT+13.)*TPl/24. 407' 

C ( 4 ) = T 2M*TP1 / 16* 408 

15 L = L-l 409 

16 M = X I NTF ( YC ( 5 ) ) - 2 410 

N = XINTF(YC(6) ) - 2 411 

XX = ( - I NTF { YC ( 5) )+YC( 5) )*2.0 - 1.0 412 

YY = ( — INTF(YC(6> ) +YC ( 6 ) ) * 2 • 0 - 1.0 413 

XP1= XX + 1.0 414 

XM1 = XX - 1.0 415 

YP 1 = YY + 1.0 . 416 

YM 1 = YY - 1.0 417 

X2M= XP 1*XM 1 418 

Y2M= YP1*YM1 419 

P < 1 ) = -XM1*X2M/16. 420' 

P ( 2 ) = ( ( 3.*XX+2. ) *XX-9. )*XM1/16. 421 

P ( 3 ) = ( -XX*XX+9 • ) / 8 • - P ( 2 ) 422 

P ( 4 ) = XP1*X2M/16 . 423 

0(1) = -YM1*Y2M/16. 424 

Q ( 2 ) = ( ( 3.*YY+2. ) *YY-9. ) *YM1/16. 425 

Q ( 3 ) = (-YY*YY+9. ) /8. - 0(2) 426 

0(4} = YP1*Y2M/16. 427 

IF (LR) 17*18*17 ’ / 428 

17 PX ( 4 ) = ( 3.*XX-1. )*XPl/8. 429- 

PX( 1 )= XX/2. - PX < 4 ) 430 

PX ( 2 ) = (9.*XX-ll.)*XPl/8. 431 
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PX ( 3 ) = -XX/ 2 • “ PX ( 2 ) 432 

QY ( 4 ) = (3.*YY“1* )*YPl/8* 433 

QY ( 1 ) * YY/2. - QY ( 4 ) . 434 

QY ( 2 ) = (9.*YY“11. )*YPl/8. 435 

QY t 3 ) = -YY/2* - QY ( 2 ) • 436 

18 DO 27 K= 1 » K4 437 

HP ( K ) = 0.0 438 

CP ( K ) = 0.0 439 

SP ( K ) = 0.0 440 

IF (LR) 19.20.19 441 

19 HPX ( K ) = 0.0 442 

HP Y ( K ) = 0.0 443 

CPX ( K ) = 0.0 444 

CPY ( K ) = 0.0 445 

SPX ( K ) = 0.0 446 

SP Y ( K ) = 0.0 447 

20 KK = ( (K+L)*32+N)*18 + M - 594 448 

DO 23 J=1 » 4 449 

HD ( J ) * 0.0 450 

CD ( J ) = 0.0 451 

SD(J) = 0.0 452 

IF (LR) 21,22.21 453 

21 HS ( J ) * 0.0' 454 

CS(J> = 0.0 455 

SS(J) = 0.0 456 

22 JJ = 18* J + KK . 457 

DO 23 1=1,4 458 

II = I + JJ 459 

HT ( I ♦ J ) = XHD ( I I ) 460 

ID = MHD(II) .AND. MASK , 461 

ID = ID * 16777216 462 

CT ( I » J ) = COSF(ARG) r 463 

23 ST ( I . J ) = SINF(ARG) 464 

DO 25 1=1,4 . 465 

DO 25 J = 1 ♦ 4 466 

HD ( I ) = Q ( J ) *HT ( I , J ) + HD ( I ) 467 

CD ( I ) = Q ( J ) *CT ( I ♦ J ) + CD ( I ) 468 

SD ( I ) = Q ( J ) *ST ( I , J ) + SD( I ) 469 

IF (LR) 24,25,24 47C 

24 HS ( I ) = P ( J ) *HT ( J , I ) + HS ( I ) 471 

CS ( I ) = P ( J ) *CT ( j » I ) + CS(I) 472 

SS(I) = P ( J ) *ST ( J ♦ I ) + SS(I> 472 

25 CONTINUE 47 *■ 

DO 27 1 = 1,4 " 475 

HP ( K ) = HD ( I ) *P ( I ) + HP ( K ) 47( 

CP ( K ) = CD ( I ) *P ( I ) + CP ( K ) , 47* 

SP ( K ) = SD( I ) *P ( I ) + SP ( K ) / 47( 

IF ( LR ) 26 ♦ 27 , 26 47< 

26 HPX(K)=HD( I )*PX( I ) + HPX(K) 48( 

CPX (K) =CD( I )*PX( I ) + CPX(K) 48 J 
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SPX(X)*SD( I )*PX( I ) + SPX(X) 482 

HPY ( X ) *HS ( I ) *QY ( I ) + HPY(X) 483 

CPY ( X ) *CS ( I ) *QY ( I ) + CPY(X) 484 

SPY(X)«SS( I )*QY( I ) + SPY < X ) 485 

27 CONTINUE 486 

H « 0.0 487 

CX - 0.0 488 

SX = 0.0 489 

IF (LR) 28.29*28 490 

28 HX « 0.0 491 

HY « 0.0 492 

CXX « 0.0 493 

CXY * 0.0 494 

SXX * 0.0 495 

SXY « 0.0 496 

29 DO 31 X=1 * X4 497 

H = C ( X ) *HP ( X ) + H 498 

CX = C(X)*CP(X> + CX 499 

SX = C(X)*SP(X> + SX 500 

IF (LR) 30,31*30 501 

30 HX « C(X)*HPX(X) + HX 502 

HY = C ( X ) *HPY ( X ) + HY 503 

CXX « C(X)*CPX(X) + CXX 504 

CXY = C ( X ) *CPY ( X ) + CXY 505 

SXX = C(X)*SPX(X) + SXX 506 

SXY = C( X ) *SPY ( X) + SXY 507 

31 CONTINUE 508 

RAD » SORTF ( CX*CX + SK*SX) 509 

CX = CX/RAD 510 

SX = SX/RAD 511 

RETURN 512 

END 513 

SUBROUTINE ANGLE 514 

DIMENSION MHD ( 4608 ) * YC ( 8 ) *YVARS( 8 ) 515 

COMMON MHD* T C*YC*H*HX»HY *CX*CXX»CKY *SX*SXX*SXY »XX*YVARS*XLG*XLT * 516 

+ A * B* CC »DA»DB * DC »LR . _ ‘ * 517 

DELX = YVARS ( 5 ) - 24. 518 

DELY = YVARS ( 6 ) - 16. 519 

COSXX= -DELX*CX - DELY*SX , 520 

SINXX 3 DELX*SX - DELY*CX 521 

IF (COSXX) 2.1*2 522 

1 XX = SIGNF(90.»SINXX) . . 523 

GO TO 6 524 

2 XX = ATANF(S INXK/COSXK. )*57 *.29377951 525 

IF (COSXX) 3.6*6 _ 526 

3 IF (SINXX) 5*4,4 527 

4 XX = XX + 180. ' 528 

GO TO 6 7 • 529 

5 XX - XX - 180. 530 

6 IF (XX) 7*8*8 531 
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7 XK = 360. + XK 532 

8 XT * DELX*. 98400775 - DELY*. 17364818 533 

YT « DELX*. 17364018 + DELY*.9B480775 534 

RAD = SQRTF(XT*XT + YT*YT ) 535 

IF (XT) 10,9.10 536 

9 XLG = SIGNF (90*0. YT ) 537 

GO TO 14 538 

10 XLG « ATANF(YT/XT >*57.29577951 539 

IF (XT) 11,14,14 540 

11 IF (YT) 13.12.12 541 

12 XLG * XLG + 180. 542 

GO TO 14 543 

13 XLG « XLG - 100. 544 

14 XLT * — ATANF ( RAD/ 31. 205 )* 114. 591559 + 90.0 545 

RETURN 546 

END 547 

SUBROUTINE AP3 548 

DIMENSION MHD ( 4608 ) , YC ( 8 ) , YVARS( 8 ) 549 

COMMON MHD.TC»YC,H,HX.HY»CK,CKX»CKY*SK.SKX*SKY*XK*YVARS,XLG*XLT* 550 

+ A » B » CC » DA ♦ DB . DC . LR 551 

R1 = SQRTF( ( .062760850324*H-.60018313990)*H+4. 7014047597) 552 

VF = 0.02 154199761 9* H + 19.270272298 - R1 553 

R2 = SQRTFU .060104035000*H-.96636105838 )*H+6. 1294779871) 554 

B = -0.12663045716*H + 19.505778258 - R2 555 

IF (H-17.) 1,1,2 556 

1 R3 =SQRTF( (• 083601632403*H-1. 3340008783 )*H+7.1705253492 ) 557 

VH = — 0 . 2479 165 049 0*H + 19.793624009 - R3 558 

GO TO 3 559 

2 R4 =SQRTF ( ( • 05 5 777 5 3 32 14*H- 3. 0851911409 ) *H+45# 69817 0763 ) 560 

VH * -0.31013284640*M + 14.848653764 + R4 561 

3 A = ( VF+VH ) * .5 562 

CC = A-VH 563 

IF (LR) 4.8.4 564 

4 DVF= ( - .06276 0050324 *H+. 30009156995 ) /R1 +0.021541997619 565 

DB = (-.0 60 1 04035000+H+ .48 318052919 )/R2 - 0.12663045716 566 

IF (H-17.) 5.5.6 * 567 

5 DVH C ( -• 0836 01 63240 3*H+. 6670004391 5 ) /R3-0 • 24791650490 568 

GO TO 7 569 

6 DVH= ( .055777533214*H-1.54259557045)/R4-0. 31013284648 570 

7 DA =(DVF+DVH)*.5 571 

DC =DA-DVH ' 572 

8 RETURN 573 

END 574 

END 575 

FINIS 576 

EXECUTE. 

OB 0574 BLEICK BOX 6 

C2AP3 DECEMBER 6 1966 JUL26,66 1 

666JUL26. 66 J 207 00.-122.5 37.9 139.6 35.6 10024 7 
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SUBROUTINE AP2 

DIMENSION MHO (4608 ) *YC(8) »YVARS(8) 

COMMON MHD.TC*YC*H.HX*HY»CK*CKX*CKY*SK*SKX*SKY*XK*YVARS*XLG*XLT* 

+ A * B * CC »DA *DB » DC »LR 

R1 = SQRTF ( ( .041783709356*H-. 42321401072 ) *H+2. 2342337579 ) 

VF » -.028281950577*H + 17.494735347 - R1 

R2 » SQRTF ( ( ♦ 05 8458667266*H-. 90 72944906 5 ) *H+6. 4033842487 ) 

B “ -0 * 1452000 1518 *H + 18.530490911 - R2 
IF (H-15.) 1*1*2 

1 R3 “SQRTF ( ( .23341292994*H-3.1096617758)*H+29. 275404601) 

VH » -0.2515261 43 53*H + 21.408836894 - R3 

GO TO 3 

2 R4 “SQRTF ( ( • 14668786198*H-6. 8828319323 ) *H+105. 12448592 ) 

VH “ -0. 369702342 18*H + 11.346369501 + R4 

3 A “ ( VF+VH ) * • 5 
CC “ A-VH 

IF (LR) 4*8*4 

4 DVF“ ( -.04 17837093 56*H+. 21 160700536 ) /R1 - .028281950577 

DB » ( -.058 4586672 66*H+. 45364724533 ) /R2 - 0.14520001518 
IF (H-15.) 5*5*6 

5 DVH= (—.23341292994*H+ 1.554830 8879 )/R3- 0.2515261 43 53 
GO TO 7 

6 DVH“ ( .14668786198*H-3. 44141596615 )/R4-0. 36970234218 

7 DA “ (DVF+DVH)*.5 
DC “DA-DVH 

8 RETURN 
END 

—COOP » BOX 6»BLEICK* I /1/0/2/S/ 56/57 *3*1 00 00*0* TAPE - 6 DEC 66. 

-FTN »L *A » E. 

PROGRAM TAPE 0 

DIMENSION X(900)*Y(900) *MD (63*63 ) *ND ( 3969 ) *MHD ( 18*32*8) 1 

EQUIVALENCE (MD*ND) * ( ID»ARG) * (TEMP* ITEM) ♦ ( IH*H) 2 

REWIND 1 .3 

REWIND 2 4 

ISCALE = 2 00000 000 0000 000 B ' 5 

MASK = 77 7777 770000 00008 6 

C READ COORDINATES FOR MAP GRID OF DRAW SUBROUTINE 

READ (50*1) ( X ( I)»I“1*390), ( Y(I ) * I “1 *390 ) 7 

1 FORMAT (15F5.3) 0 

WR I TE ( 2 ) X * Y 9 

IF ( I OCHECK * 2 ) 2*4 10 

2 WR I TE ( 51 *3 ) 11 

3 FORMAT ( 37H0 PARITY ERROR OCCURRED ON X*Y WRITE/) 12 

4 DO 43 K= 1 ♦ 8 13 

C READ WAVE DIRECTION FROM FLEET NUMERICAL WEATHER FACILITY TAPE 

BUFFER IN ( 1 *2 ) ( ND ( 1 > *ND( 3969 )) 14 

5 IF (UNI T * 1 ) 6*14.8*10 7 15 

6 GO TO 5 16 

8 WRI TE ( 51 *9 ) K 17 
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c 



c 

c 



c 



9 FORMAT (44H0 DIRECTION EOF OR EOT ERROR OCCURRED ON K=I1/) 18 

STOP 19 

10 WR I T E ( 5 1 * 1 1 ) K 20 

11 FORMAT ( 4 9 H 0 DIRECTION PARITY OR LENGTH ERROR OCCURRED ON K*I1/) 21 

M = LENGTHF ( 1 ) 22 

IF (M-3969) 12*14,12 23 

12 WRITE(51»13) M 24 

13 FORMAT ( 2 8 H 0 DIRECTION BUFFER LENGTH =15/) 25 

COMPUTE WAVE DIRECTION K FROM X AXIS OF STEREOGRAPHIC GRID 



14 DO 42 1=1,18 

DELX = 1-24 
DO 42 J= 1 » 32 
DELY = J-16 

ID = MD( 1 + 8, J-t-16) /2048 + ISCALE 
ARG = (ARG + 0.0)*11. 17010721 
COS = COSF(ARG) 

SIN = SINF(ARG) 

ABS = -DELX*COS - DELY*SIN 
ORD = DELX#SIN - DELY*COS 
IF (ABS) 36,35,36 

35 TEMP = SIGNF( 1.5707963260, ORD) 

GO TO 40 

36 TEMP = AT ANF ( ORD/ABS ) 

IF (ABS) 37,40,40 

37 IF (ORD) 39,38,30 

38 TEMP = TEMP + 3.1415926536 
GO TO 40 

39 TEMP = TEMP - 3.1415926536 

40 IF (TEMP) 41,42,42 

41 TEMP = TEMP + 6.2831853072 



26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 



SHIFT FLOATING-POINT K TO LOWER HALF OF MHD MATRIX WORD 

42 MHD ( I » J , K ) = ITEM/16777216 47 

READ WAVE PERIOD FROM FLEET NUMERICAL WEATHER FACILITY TAPE (NOT USED) 
BUFFER I N ( 1 , 2 ) ( ND ( 1 ) , ND( 3969 ) ) 48 

15 IF (UNIT, 1) 16,24,10,20 , 49 

16 GO TO 15 50 

18 WR I T E ( 5 1 , 19 ) K 51 

19 FORMAT ( 4 1H0 PERIOD EOF OR EOT ERROR OCCURRED ON K=I1/) 52 

STOP 53 



20 WRITE (51,21) K 

21 FORMAT ( 46H0 PERIOD PARITY OR LENGTH ERROR OCCURRED ON K=I1/) 
M = LENGTHF ( 1 ) 

IF (M-3969) 22,24,22 

22 WR I TE ( 5 1 , 23 ) M 

23 FORMAT (25H0 PERIOD BUFFER LENGTH =15/) 

READ WAVE HEIGHT H FROM FLEET NUMERICAL WEATHER FACILITY TAPE 

24 BUFFER IN(1,2) ( ND ( 1 ) *ND( 3969 ) ) / • 

25 I F ( UNI T ♦ 1 ) 26,34,20,30 

26 GO TO 25 

28 WR I TE ( 5 1 , 29 ) K 



54 

53 

56 

57 

58 

59 

60 
61 
62 
63 
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29 

30 

31 



32 

33 



FORMAT (41H0 HEIGHT EOF OR EOT ERROR OCCURRED ON K =» 1 1 / ) 

STOP 

WR ITE(51*31) K 

FORMAT (46H0 HEIGHT PARITY OR LENGTH ERROR OCCURRED ON K=I1/) 
M = LENGTHF ( 1 ) 

IF CM-3969) 32*34* 32 
WR I TE ( 5 1 ♦ 33 ) M 
FORMAT (23H0 HEIGHT 



BUFFER 

PACK FLOATING-POINT H AND K IN 



LENGTH * 15/) 
MHD MATRIX 



34 



43 



44 

45 

46 



1=1*18 
J = 1 * 32 

MD( 1 + 8 * J+16 ) /2048 
(H + 0.0)*64. 

IH .AND. MASK 

= MHD( I »J»K ) 



44*46 



+ iscale 



+ IH 



DO 43 
DO 43 
IH = 

H = 

IH = 

MHD ( I ,J,K) 

REWIND 1 
WR I TE ( 2 ) MHD 
I F ( I OCHECK ♦ 2 ) 

WRITE( 51*45) 

FORMAT (37H0 
END FILE 2 
REWIND 2 

PRINT PART OF MHD' MATRIX TO CHECK PACKING OF DATA 

WRITE! 51*47 ) ( ( ( MHD( I*J»K)»I=1*7)*J=1*32.31)»K=1»8) 
47 FORMAT (7017/) 

STOP 

END 

END 

FINIS 

:xecute. 



PARITY ERROR OCCURRED ON MHD WRITE/) 



64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 



978 


978 


17461 


17461 


978 


978 


2070 


4931 


4133 


3388 


2698 


2065 


1491 


978 


978ABS 


1 


491 


2064 


2697 


5198 


4550 


3956 


3418 


2937 


2515 


2152 


1851 


1612 


1435 


1321 


1271ABS 


2 


284 


1361 


1501 


1704 


1969 


2296 


2684 


3130 


3635 


4197 


4814 


5484 


6205 


6976 


7794ABS 


3 


764 


9841 


8966 


8144 


7378 


6669 


6021 


5437 


4917 


4465 


4082 


3769 


3528 


3359 


3262ABS 


4 


240 


3290 


3414 


3611 


3879 


4219 


4628 


5106 


5651 


6259 


6931 


7661 


10146 


9307 


8535ABS 


5 


834 


7207 


6658 


6188 


5800 


5496 


5278 


5146 


5102 


5145 


5275 


5491 


5793 


6179 


6647ABS 


6 


196 


7821 


8520 


9291 


101281 


1028 


119871 


3000 


14062 


17461 


17461 


1647715523 


14600 


13714ABS 


7 


8681 


2067 


11312 


10608 


9957 


9363 


8827 


8353 


7942 


7596 


7316 


7104 


6960 


6885 


6880ABS 


8 


944 


7078 


7280 


7550 


7886 


8288 


8753 


9279 


98651050711204 


1195112746156721482RABS 


9 


0251 


3265 


12551 


11887 


112771072310227 


9793 


9422 


9116 


8877 


8705 


8601 


8567 


8601ABS10 


, 70 5 


8877 


9116 


9422 


9793102271072311277118871255113265 


140261482915673 


16551ABS11 



<f6 11 7461 16663158971 516714475138251 3220 12664121 591 1708 1131 11097 3 106941 0475 4BS12 
3181022410193102241031810475106941097311311 11708 12 159 1266413220 1382 514475ABS13 
1671 5897166631746 11746 11682716221 15645 151 00 14590141 17136821 3288 1293612628ABS14 
.361 2 1491 1979 1 1 8581 1785 1 1 7601 178 5 118 58 11 979 12 149 12 366 12 628 129 361 328 8 13682A5S1 5 
1171 4590 15100 15645 1622 11 682717461 17461 1698916539 161 13 157 121538 8 1499 114674ABS16 
18714131 13908 137171356 01 3438 133501 3297132791 32971 3350 13438 13560 1371 71 3908ABS17 
31143871467414991153381571216113165391698917461 17461 17 150 168 551 657 71 631 7ABS18 
*741585 11 5647 15462 152991 5 156 1503 51 4935 148 57 14801 14768 1475714768 1480 11 43 57ABS19 
351 5035 15 156 152991 5462 15647 15851 16074163 17 16577 168 55 171 50 17461 1746 11 73 12 ABS20 
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1 71 72 1704 1 169 19 16806 16702 16600 1652 3 16440 16 382 16327 1628 116 246 16220 1620 5 16200ABS;i 
16205 16220 162461628 11632 7 16382 16448 16523 16600 16702 16806169 19 1704 11 7172173 12ABS| 
17461174611672215277174611746113665110241746117461 9667 70621746117461 3797ABS3 
978 9781746117461 978 9 78 1 746 1 1 746 1 978 9781746117461 978 97 8 1 746 1 ABSift 

17461 978 9781746117461 978 97 017461 17461 978 9 781746117461 2520 5878ABS;5 

1746117461 86421099117461174611304016870 1746 11746116537 ABSiS 

6292870528705 629 629 2224 629 629 1612 2637 3700 4797 5927 70 8 52 53 2 20R D 1 

264802 760928 70528 7052772 126 703256552457823478 22 3572 12 17200 63 18898 1772 5 165470RD? 
1536814192130221186110713 9500 0467 7377 6312 5275 4271 3302 2369 1478 6290RD3 
629 1408 2 2- '+2 3127 4061 5039 6059 7116 8207 9328 10473 1 1640128241 40201 52240RD ft 
16432 1763918841 2003321 2 112 2371 23507246 1725695 2673927743287052870 52777 2267830RD| 
2574224655235282236421171199531071017470162161496213714124781126110067 89030RD5 
7774 6686 5645 4654 3720 2046 2037 1297 629 629 976 1379 1844 2372 29580RD7 

3601 4298 5046 5842 6602 7562 8400 943 1 1041 0 1 1 41 5 1 2 440 1 348 1 14 5 33 1 5593 166 560R D 3 

1771618770 19 81 320 841 21 8492 2 032 2 3708 2 47102 5 597 26443 2 72 462 8001 287052 870 52 8 1630RD9 
275642691 126207254562466023025229532205021 11 82 0164 191 90 18203 172 05 162031 520 20R DO 
1420413217122431128810357 9453 8502 7746 6951 6200 5496 4843 4244 3702 32190RD1 
2797 4582 5005 5483 6013 6593 7220 7890 8602 9350 10 1 3 2 10943 1 1 78 0 1 263 8 1 3 5 1 40R D ? 
1440315301 16204 17 106 18004 10093 19769206272 14642 2 2762 305723 806245 17 251 8 72 58 140RD 3 
2 6 3942 69242 7402 2 78 2 526052 2 5 6602 52 272 475624248 2370 52 3 130 22 5 2 52 189 3 2 123 7205 590RD4 
1986219149184241768916948162041545914710139831325812545118491117010514 98820RDE 
9277 8702 8159 7651 7100 6747 6355 0154 8507 8807 9294 97261018 1 106571 1 1 540RD6 
11 668 121 99 12 744 1330 113869 14445 1502 8 156 15 162041 6792 17379 17962 185 38 191 061 96630RD7 
20 2082 07392 12 542 17 50222262268 123 11 32 3520 23900 24253 2 2 360 220 602 1743 2 1412210670R.De 
20 70920 33 8 19957 19 566 19166 187 57 10 342 1792 11 7496 17067 16636 162 041 577 11 53401 49 110RDS 
14486 1 4065 1 3650 13242 1284 11 245012 069 11690 11 340 10995 1066410347 10047 1222 2 124420R DC 
1266712897 13 133 13372 1361 61 3064141 151437014627 14086 15147 154 10 156741 593 8 162040RD1 
16 469167 3316997172601752117781180301029210 54310 791 1903 5 192 7 5195 1 0 197401 99660RD- 2 
20 1852 726 128705287052 52242 360520 7052 8705 224662 1462 2870 528 70520608 198642 870 50RD? 
28 7052792 9 19200 185942 5558 2 3 3571 803 11 75002 1277 19278 16989 16491 1732 71 53941 59970RDi 
1549913448114591499114464 9395 72161390713309 4876 23181265511926 629 6290RD5 

1109310117 629 629 8930 7461 629 629 5521 2814 629 0RD6 



List of cards in programs TAPE, VC2AP3 and AP2 requiring 
changes if the geometric dimensions and/or origin of the MHD 
array are changed: 

TAPE : 1, 26, 27, 28, 29, 30, 72, 73, 74 and perhaps 85, 

and the 52 input BCD data cards for the computer 
plotting of the map grid. 



VC2AP3 : 


1 , 


86, 87 


, 88, 


89, 




157 


, 158, 


163 , 


1 64. 




303 


, 304, 


33 1 , 


332, 




519 


and 5‘ 


49. 




AP2 : 


DIMENSION 


statement 



05, 106, 109, 110, 


141 , 142, 


202, 


203 , 


223 , 


224, 


296, 


297, 


346 , 


348, 


448, 


457, 


51 5, 


518, 
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IDENTIFICATION: 

TITLE: General Graph Output Subroutine 



CO-OP ID: 

CATEGORY: 

PROGRAMMER: 

DATE: 

REVISED: 



J7-NPS-DRAW 

Output for Off-Line Plotting 
J . R. Ward 
February 1964; 

June 1965 



(FORTRAN 60) 



B. PURPOSE: 



This subroutine, when provided with the necessary information, 
generates a magnetic tape in the proper format for subsequent 
off-line graph plotting using the CDC 160 or 160A GRAFPLOT program 
and a CalComp 165 Plotter (see references 1 and 2). Provision is 
made for curve drawing and point plotting, automatic scaling, graph 
titling and axis annotating. An attempt was made to provide a 
considerable amount of flexibility, at the expense, necessarily, . of 
a relatively large number of arguments and a rather high memory 
requirement. ■; > . - * »•.“ i 



C. USAGE: 

1. Definitions: 



In what follows the word "graph" will be taken to mean one 
piece or frame of graph paper on which there may be plotted 
one or more curves and/or sets of points. A "curve" will 
mean a continuous line generated by the sequence of straight 
lines joining successive points of the set defining the curve. 

A "point plot" will describe the representation of a succession 
of points by means of symbols (such as a cross) on the graph. 
The points are not connected in a point plot. 



2. Calling Arguments: 



All necessary information is transferred to DRAW through the 
calling arguments. The call statement is: CALL DRAW (NUMPTS, 

X , Y , MODCURV , ITYPE, LABEL , ITITLE , EXSCALE , YSCALE , IXUP , IYRIGHT , 
MODEXAX , MODEYAX , IWIDE , IHIGH , IGRID , LAST) 

It is important to realize that one and only one curve or set 
of points is plotted each time DRAW is called. However, it is 
possible to call DRAW repeatedly if several curves and/or sets 
of points are wanted on one graph. 






/ . 



dr' 



.J j<' 
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/ 



I 



The calling argupjontu avo as follows: 

■ ■ f ; 

a. NUMPTSJ Tho number of points defining a curve (2 ‘s NUMPTS 

> 4 . 

)'< £ 900) , or the number of points to be point plotted 
* (2 NUMPTfl £ 30<). 



b. 



c. 



d. 



t 






Os 



X: 



Y: 



MODCURVj 

■ 0 
- 1 

- 2 

- 3 

ITYPE: 

■ 0 

■’ 1 

■ 2 

- 3 
■ A 

- 5 



t The array of X-ordinates < | i 1<V^ for 1-1,2,..., 

NUMPTS) . X must be dimensioned at least equal to 

< > ' 

"" humors in the calling program. 

All points will be considered to have the same 

X-ordinate if |X - X JilO' '. The common value 
*ax min 

will be put equal to aero if IX I £ 10 
Tho array of Y-ordinates, with properties corresponding 
to the X-ordinates, above. Y must be dimensioned at 
least equal to NUMPTS in the calling program. 

J x 

Controls the number of curves, and/or sets of points 
■ v on one graph* 

\ 

This is the only curve, or set of points, to be* plotted 
on this graph. 

This is the first of two or more curves, and/or sets 

of points, to be plotted on this graph. 

This is an intermediate curve, or set of points. 

» » 

This is the last curvo, or sot of points, for this 

» 

graph. 

i 

Controls the type of plot (l.e., curve or point plot): 

* t- 

This set of points is to be represented by a curve . 

i 

Those points are to be plotted with a cross (x). 

These points are to be plotted with a plus (+). 

These points are to be plotted with a square (□). 

These points are to be plotted with a diamond ( 0 )* 

These points are to be plotted with a triangle (A). 



This is a Hollerith curve or point identifier. If a 
curve is being drawn, LABEL must have 4 characters 
(including any blanks), and these will be reproduced 
beside the end of the curve. This argument can be 

f 

set in the calling program by a statement such as 
LABEL - 4 H a 0NE - • (^“ blank) 

or LABEL = 4H1234 

or LABEL ■ 41k. AAA The latter must be used 

when no label is wanted. 

If a set of points is being plotted, LABEL is an 8-character 
identifier. The first 4 characters will be reproduced 
beside the first point, and the last four characters 
will be reproduced alongside the last point. This 

argument can be set by statements such as 

. 4 . ... 

LABEL = 8HFRSTLAST 
or LABEL - 8H AAAAA 0NE 

or LABEL = 8 H a ONE a 123 

or LABEL *> 8 Haaaaaaaa The latter must be used 

if no labels are wanted. 

The above arguments, a. through f. (and q.), have 
meaning every time DRAW is called. On the other hand, 
the remaining arguments, g. through p., have no meaning 
except when MODCURV ■ 0 or l. v ' r- J 



g. ITITLEJ An array of twelve 8-character Hollerith words, the 

first six of which will form the first title line, 
and the last six the second. The array must be dimen- 

• i 

sioned 12 in the calling program , must contain the 

. • '. 4 .. 

user's job identification, and must have unwanted 
• characters set to blank. For example: 

DO 1 I - 1,12 -i ' •'•>■■■' ■ 

x ’ » ** ■' 

- V 1.; 'ITITLE(I) - 8H ■’ v ' ' 

ITlTLE(l) - 8 H a SMITH, a 
ITITLE(2) - 8 HJ. a J.aaa 

ITITLE<7) - 8 H a TESTIT. 

_qq 99 

h. EXSCALE: X-scale in units per inch (10 a . EXSCALE 10 • ) . 

EXSCALE will always be rounded off to one figure 
significance. If EXSCALE ■ 0, the X-scale will be 

i • . •> 

computed by DRAW. This is called auto-scale. . 

i. YSCALE: Y-scale in units per inch, with properties correspond- 

ing to those of EXSCALE. 

. * 

.* < \ 

jv IXUP: < Distance, in inches, of the X-axis from the bottom 

of the graph (0 £ IXUP 6 . IHIGH) . This argument will 
be ignored unless MODEXAX =2, see below.. 

k. IYRIGHT: Distance, in inches, of the Y-axis from the left of 

the graph (0 — IYRIGHT £ IWIDE) . This will be ignored 

t 

unless MODEYAX ■ 2, see below. / 

l. MODEXAX: Determines the mode of the X-axis location: 

■ 0 The X-axis will be located automatically by DRAW, 

with the origin of Y on the graph. 



- 1 



• • -i '> 

- 2 

m. HODEYAX: 

i • 

n* I WIDE: ^ 

4 r.v. *: 

o»- 

o. IHIGH: 

' ;;; r,v 

p. IGRID: 



LAST: 



.,‘Z s • a 

. / r r* ; • * j 

« v : jt i j 

♦ • i . . • M • , , 

p Jf 

< Ih ■’/ 

, »; J) i, i > „■ 



• . . 

. '■ r'Vo 



\. { O 



.r-< 



»}■.*. .. 
7 ■» 



• The X-axis will be automatically located by DRAW, 
with the origin of Y moved (in one's imagination) an 
integer number of inches above or below the graph, 
if this is appropriate. This option can be used only 
if the Y-scaling is automatic (YSCALE ■ 0) . 

i 

The X-axis location will be as specified by IXUP. 
Determines the mode of Y-axis location in the same 
way as MODEXAX, above, governs the X-axis location. 

Width of graph in inches (1 £ IWIDE £9). If IWIDE 
is out of this range, a value of 8 will be assumed. 

Height of graph in inches (1 ^ IHIGH 6: 15) . If IHIGH 
is out of range, a value of 8 will be assumed. 

If IGRID “ 1, a 1" x 1" grid will be superimposed on 
the graph. This is useful only if plain paper is used 
on the CalComp Plotter. 

Indicates to the calling program whether the previous 
plot was completed successfully. The codes are: • 

0 Last plot was completed successfully 

'* 

■ 1 Last plot was not completed successfully. 

y ... 

■ 2 Last plot was not completed successfully, and 

no further graph output will be . attempted until 
DRAW is next entered with MODCURV «* 1 or 0. 

Si ■ 3 An attempt was made to enter DRAW with MODCURV 
[ > -i' M or 0 while the error lockout was set. 

■ * • t>\, ■ *, 

A \ ' 

This argument must always be a variable, name and 



• • •' 

. , tl, ,• 

j;,\! ? 



... i * 



never a number in the call statement. 
-41-^ • ‘ ' 



i 



3* Notes and Comments: 

a. The graph scales and, if MODEXAX - 1 and/or MODEYAX - 1, the 
amounts of origin offset are always output as part of the graph title. 

b. Each time a graph is completed, a message to this effect is printed 
on both the standard output and the console typewriter. 

c. There are internal checks of the input to DRAW to prevent incorrect 
use. If an input error is detected, an attempt will be made, where 
possible, to complete the plot. If an argument is "corrected" in 
this process, the user will be so informed on the standard output. 

If it is not possible to complete the plot, the user will be informed 
of the reason by a message on the standard output. 

d. If part or all of a curve would fall more than 0.6" laterally 
beyond the ends of the X-axis, or 0.7" vertically beyond the ends of 
the Y-axis, the X and/or Y ordinates will be limited so that the curve 

. will typically become a line along part or all of the boundary of 
the graph as here defined. 

e. If one or more points of a point plot would fall outside the graph 
area, the plot of that point, or points will be inhibited. The 
number of such points will be reported to the user on the standard 
output. 

f. It should be pointed out that the X and Y scaling and axis locating 
processes are entirely independent, so that, for example, X might 

be auto-scaled, while the Y-scale is specified. At the same time 
the X-axis might be located automatically, while the Y-axis location 
is specified. 

g. It must be remembered that the scales and axis locations of a 
multi-plot graph are set when DRAW is called for the first time 
(with MODCURV ■ 1). Thus the user must attempt, at that time, 
to achieve scaling and axis location which will be appropriate 

to all the plots he intends to make on the one graph. Particularly 
if the automatic features of DRAW are selected, foresight will be 
demanded of the user in this respect. 

4. Auto-Scale Properties: 

The scale factor is chosen from amongst the values 1,2 or 5 units per 
inch, or some power of 10 times one of those values. A curve, or 
set of points which is plotted with auto-scale will normally lie entirely 
within the graph area as defined in 3.d., above. The only exception may 
occur if an axis is placed, by the user, along one edge of the graph 
(e.g., IXUP *» 0, MODEXAX «* 2). In such a case, points "outside" 
the axis are not considered in the selection of a scale factor (e.g.,’ 
negative X^ do not affect the choice of scale when IXUP ° 0) . If 
automatic axis location as well as auto-scale is selected, the plot, 
if it does not fill the graph area, will be placed as far as possible 
towards the bottom-left of the graph area consistent with the fact that 
the axes can be set in increments of 1" only. 

-42- 



5. Space Required: 3960 cells (excluding the input arrays). 

f 

6. Temporary Storage: None 

7. Error Print-Outs: There are a large number of possible error print-outs. 

These are all self-explanatory. 

8. Error Returns: All error returns are preceded by self-explanatory error 

printouts. An error indication is transferred back to the calling 
program through the argument LAST. 

9. Error Stops: None. 

'10. Tape Mountings: Logical Tape //8 will receive the binary graph output. 

The standard monitor output will receive the messages to the user. 

11. Output Format: The format of the binary graph output records on 

magnetic tape is described in reference 1. The only difference is that 
in this program the interpolation option is by-passed (set to zero in 
the graph output record)’. See reference 2. 

12. Selective Jump arid Stop settings: None. 

13. Timing: Variable, depending upon the number of points and the options 
chosen. Typically less than one second per curve or point plot. 

14. Accuracy: The accuracy of results is equal to the resolution of the 

CalComp Plotter, that is, 0.01" in both the X and Y directions. 

15. Equipment configuration: CDC 1604 with FORTRAN 60 compiler and 

Library. A CDC 160 or 160A with CalComp 165 Plotter is needed for the . 
off-line plotting using the appropriate GRAFPLOT program. 

fi. REFERENCES: 

1. Weir, Maurice D. , Spritzer, Milton and Mclihenny, D. W. , "160-A 
Graph Plot Program," Ident *2001, SWAP Library, 15 August 1962. 

2. Hogg, R. L. and Glover, D. C. , "160 Grafplot Routine," Writeup 
available from Computer Facility, U. S. Naval Postgraduate School, 

1. April 1963. 

3. U. S. Naval Postgraduate School, Thesis, "Control System Programming, 
Remote Computing and Data Display," by Robert Lee Hogg and Dennis 

C. Glover, 1963. 



1J.B,. References 1 and 2 are included in Reference 3 as Enclosures 
2 and 1, respectively. 
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